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Dissolved oxygen in sea water is a major factor affecting marine habitats and bio¬ 
geochemical cycled^®! Oceanic zones with oxygen deficits represent significant por¬ 
tions of the area and volume of the ocean^and are thought to be expandin^H^ The 
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Peruvian oxygen minimum zone is one of the most pronounced and lies in a region 
of strong mesoscale activity in the form of vortices and frontal regions, whose effect 
in the dynamics of the oxygen minimum zone is largely unknown. Here, we study 
this issue from a modeling approach and a Lagrangian point of view, using a cou¬ 
pled physical-biogeochemical simulation of the Peruvian oxygen minimum zone 
and finite-size Lyapunov exponent fields to understand the link between mesoscale 
dynamics and oxygen variations. Our results show that, at depths between 380 
and 600 meters, mesoscale structures have a relevant dual role. First, their mean 
positions and paths delimit and maintain the oxygen minimum zone boundaries. 
Second, their high frequency fluctuations entrain oxygen across these boundaries 
as eddy fluxes that point towards the interior of the oxygen minimum zone and 
are one order of magnitude larger than mean fluxes. We conclude that these eddy 
fluxes contribute to the ventilation of the oxygen minimum zone. 

Regions of the ocean with strong O 2 deficiency in the water column are called 
Oxygen Minimum Zones (OMZs). The OMZs are differentiated in the vertical by three 
distinct layers: the oxycline (upper O 2 gradient), the core (typically with O 2 < 20 fiM) 
and the lower O 2 gradient. The Eastern Tropical South Pacific (ETSP) contains one 
of the three major permanent OMZ, with an oxycline extending from the upper 10 to 
170 meters, a core with a thickness of around 400 meterP and lower oxygen gradient 
extending to about 3700 meterP. This OMZ is maintained by the combination of signif- 
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icant rates of biological production and decomposition of sinking organic materiaP at 
the Peruvian upwelling region, and weak circulation in the shadow zone of the southern 
Pacific subtropical gyre. The circulation is then dominated by the equatorial and eastern 
boundary current system^^^. Energetic vortices, called mesoscale eddies, and filaments 
are ubiquitous in this areapl. In contrast to the oligotrophic regions of the ocean where 
mesoscale eddies can sustain biological productivitjff^, in upwelling regions stirring by 
eddies - the process of tracer gradient intensification by advection - tends to inhibit 
biological productiorU^H^. 

In the ETSP, the role played by mesoscale structures in the distribution of O 2 
within the OMZ remains unclear and we approach this issue by analyzing data from 
a coupled physical-biogeochemical high-resolution modeP^ of the regional ETSP (see 
Methods), and characterizing mesoscale transport and stirring by means of Finite-size 
Lyapunov exponent (FSLE) field^i^Ill (see Methods and Supplementary Information). 
Maxima in these fields form thin filamentary structures, the so-called Lagrangian Co¬ 
herent Structures (LCS identifying the most intense mesoscale regions and acting 
as barriers for fluid transport across them. 

In this work we focus on the transport aspects of the mesoscale-OMZ interaction, 
particularly in the OMZ boundaries, and the fluxes across them. We do not address 
specifically the biogeochemical processes occurring inside the zone which are certainly 
determinant (and are included in our regional simulation model) but we gauge instead 


3 






the physical effects of the mesoscale structures on the OMZ dynamics. This is done 
by: a) computing correlations between the (temporally averaged) O 2 concentration and 
FSLE at layers located at different depths; b) studying events of O 2 rich-waters entrain¬ 
ment into the OMZ; c) calculating the temporal average of O 2 normal fluxes across the 
northern and southern boundaries of the OMZ as a function of depth, and its correlation 
with the average mixing measurement obtained from FSLE. The outcome of these anal¬ 
yses is that, despite the important biogeochemical processes, mesoscale stirring already 
determines many important features of the oxygen distribution. 

The 20 /iM isosurface of the annual mean O 2 field for simulation year (s.y.) 21 
(Fig. [T^) gives an OMZ core with maximal horizontal extension at approximately 400 m 
depth extending between 3° S and 16° S. The higher O 2 concentrations north of 2° S are 
associated with eastward equatorial subsurface currents carrying relatively oxygen-rich 
wateili2l22l, while the southern increase of O 2 (14° S to 17° S) is adjacent to the northern 
part of the subtropical gyre. Figure [T^ also displays the annual mean backward FSLE 
field at 410 m depth, which shows a high correlation with the mean O 2 field delineating 
the limits of the OMZ core. The FSLE mean field is structured as zonal bands coinci¬ 
dent with the north and south OMZ boundaries with relatively high FSLE values when 
compared to the core region. Both bands signal stirring by the eddies formed at the 
continental shelf and advected offshore, and by other mesoscale processe^i2l22l. This 
indicates that the enhanced mesoscale activity in those areas delineates the limits of the 


4 




average OMZ core region. 


Since LCS (that we locate as maximum values of FSLE) act as transport barriers, 
large gradients of O 2 should occur across theml^. Thus we expect to find a relationship 
between the stirring intensity as measured by the FSLE and the O 2 gradient norm (see 
additional discussion in Supplementary Information; in the following, the term O 2 gra¬ 
dient refers to the norm). The relationship between both fields is quantified in Fig. 
where we plot the latitudinal variation between 18'^ S and 2° N of mean FSLE and O 2 
gradient averaged between the coast and 85°W and from 380 to 600 m depth, showing 
the coincidence in the maxima of both quantities: the maxima of ESLE indicating the 
positions of the LCS and the maxima of the O 2 gradient signalling the northern and 
southern boundaries of the OMZ. This correlation is not equally strong at all depths 
as it is shown in Eig. [T]:, where we plot the vertical profile of the Pearson correlation 
coefficient ^ , R, between the zonally averaged mean FSLE and mean O 2 gradient. 
Roughly, we can distinguish two areas in the OMZ core: a) between 190 — 350 meters 
where these quantities show correlations of alternating sign; and b) between 380 — 600 
meters where the correlation is large and positive (with an average R of 0.748). It is 
in this subregion of the OMZ core that its boundaries are strongly determined by the 
mesoscale dynamics. This is so because the OMZ dynamics is a balance between hy¬ 
drodynamic and biogeochemical processes. At these depths, the physical forcing, albeit 
lower than in the upper layers, has a variability two orders of magnitude larger than the 
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biogeochemical forcing ^ (see Supplementary Information) and its effects are clearly 
visible in the strong correlation between FSLE and O 2 gradients. 

Besides mean behaviour, individual events are also relevant, since mesoscale ed¬ 
dies are able to transport waters with different biogeochemical properties with respect to 
surrounding areas, giving rise to sporadic episodes of high O 2 patches inside the OMZ. 
Fig. 1^ shows one of these temporal sequences where an eddy dipole (with borders sig¬ 
nalled by maxima of FSLE at 410 meters) entrains water with high oxygen content (the 
red-yellow tongue at 80 — 82°W) towards the interior of the OMZ. This episode had a 
duration of approximately 3 months (9/9 to 1/12 of s.y. 21; first month displayed) and 
during this period the entrainment of these waters carried 0.4 x 10® mol of O 2 per meter 
of depth into the OMZ at this depth. Episodes of this nature are frequent at the southern 
boundary although often less intense (see Supplementary Information). At the northern 
boundary the frequency of O 2 injection episodes is higher and they last longer. The 
difference between both boundaries rests upon the spatial distribution of O 2 . Since the 
mean northern boundary is almost coincident with the large O 2 gradient region most of 
the time, any small displacement of this region will cause a significant change in the O 2 
signal right at the boundary. In the southern border gradients are more distributed and 
strong anomalous O 2 signals will be caused only by the mesoscale eddies entraining 
water across the boundary as in Fig. At both boundaries the episodic ventilation of 
the OMZ follows an offshore path consistent with the propagation of eddies and other 
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perturbations from the coastal waters to the open ocean. A characterization of FSLE 
and O 2 joint dynamics in terms of wavelet spectra, emphasizing the dominant periods, 
is presented in the Supplementary Information. 

The average amount of O 2 entering through the OMZ boundaries due to mesoscale 
processes was quantified by computing eddy fluxes of O 2 normal to the northern and 
southern limits. Small-scale turbulent diffusion produces much smaller fluxes. Eddy 
fluxes were calculated across the mean 20 /iM level boundary between 200 and 600 
meters of depth during s.y. 21, from the covariance between velocity anomalies and O 2 
concentration anomalies (see Supplementary Information). Vertical fluxes across these 
borders were always orders of magnitude smaller than horizontal ones (see Supplemen¬ 
tary Information), and thus the following discussion about horizontal components does 
also apply to the total normal flux. At the northern boundary the horizontal eddy flux 
profile is mainly positive (Fig. |^, red line), meaning that the O 2 variance due to hor¬ 
izontal eddy fluxes is bringing O 2 into the OMZ. The highest eddy fluxes are reached 
at core depths between 350 and 500 meters which is close to the depth range where the 
higher FSLE mean values at the boundary are obtained (Fig. |^, blue line), although 
the maximum value of this latter quantity (associated to the presence of the subsurface 
current^^ESl) appears deeper than the eddy flux maximum (350 vs 480 meters). Above 
300 meters the horizontal eddy fluxes are small and the minimum is obtained around 
300 meters, where the FSLE is also minimum. Globally integrated between 200-600 
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m depth and along the northern 20 jiM boundary from coast until 88‘^W, the horizontal 
eddy flow rate towards the OMZ interior is of 6.15 x Kfumol s~^, whereas the corre¬ 
sponding mean flow rate is 2 x 10^/imo/ and directed outwards (see Supplementary 
information). At the southern 20 /iM mean boundary, eddy fluxes are also positive (Fig. 
[^) along the range of depths considered, being fairly constant from 200 to 300 m, and 
nearly vanishing between 400 to 600 m depth. Integrating between 200-600 m from 
coast to 88"W, the eddy flow rate towards the OMZ interior is of 1 x nmol s~^, 
whereas the mean flow rate is 3.56 x 10^/rmo( and directed outwards (see Supple¬ 
mentary Information). 

The differences in the O 2 eddy fluxes between the northern and southern bound¬ 
aries may be understood in terms of the mesoscale activity. In the southern boundary O 2 
anomalies are caused by eddies (signalled by large FSLE) crossing the boundary. Thus 
higher O 2 eddy fluxes should be associated with higher FSLE values, which indeed is 
true looking at the profiles in Eigj^). On the northern boundary, this holds until the 
local minimum at 312 m. Below this depth, O 2 anomalies crossing the northern border 
are mainly related to, as stated above, fluctuations in the position of a large O 2 gradi¬ 
ent zone, associated to fluctuations in the subsurface equatorial currents, separating the 
subsurface O 2 rich equatorial waters from the OMZ core. 

To conclude, in this work we have addressed the role of mesoscale structures that 
populate the OMZ in the ETSP. We identified the boundaries of these mesoscale ed- 
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dies and fronts as LCS that act as barriers to transport controlling fluid interchange in 
and out the OMZ. Comparison of the FSLE approach with an exit time characterization 
(see Supplementary Information) supports this view. Despite the important biogeo¬ 
chemical processes, mesoscale stirring already shapes important features of the oxygen 
distribution. We find that mesoscale dynamics plays a dual role, which can be respec¬ 
tively associated with the average behaviour and with the turbulent fluctuations. The 
northern and southern boundaries of the OMZ core are well determined by the aver¬ 
aged mesoscale dynamics, in particular for depths between 300 and 600 meters, where 
a good correlation between mean FSLE and O 2 gradients was found. At other depths 
the relation between ESLE and O 2 may not hold, indicating significant O 2 forcing by 
biogeochemical processes. Episodic events of OMZ ventilation are produced by eddy 
stirring where waters with high O 2 content are entrained into the OMZ by the action of 
mesoscale eddies. On the whole, between 200 and 600 m depth, eddy fluxes were found 
to bring O 2 inside the OMZ at both the northern and southern frontiers, while O 2 mean 
fluxes were much smaller and in the opposite direction. The biogeochemical processes 
occurring in the interior of the OMZ would provide the dominant oxygen consumption 
sink to close the O 2 budget and maintain the OMZ core. 
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Methods 


Circulation and OMZ modeling 

The circulation and OMZ modeling in the Eastern Tropical Pacific was accomplished by 
the combination of the hydrodynamic model ROMI^ (Regional Ocean Modeling Sys¬ 
tem) and the biogeochemical model developed^ for the Eastern Boundary Upwelling 
Systems (BioEBUS). The Eastern Tropical Pacific configuration covers the region from 
4" N to 20" S and from 70" to 90" W with an horizontal resolution of 1/9" degrees (~ 12 
km) and 32 terrain-following vertical levels with variable vertical resolution (higher in 
the upper ocean). The coupled model is run in a climatological configuration previously 
validated^ for the Eastern Tropical South Pacific, and the present configuration has 
been recently validated and a sensitivity analysis performed^!. The model was forced by 
the QuickSCAT® wind stress monthly climatology and by heat and fresh water fluxes 
from the COAD^^ monthly climatology. The dynamical variables at the three open 
ocean boundaries are provided by a monthly climatology computed from the Simple 
Ocean Data Assimilation reanalysi^. Eor the biogeochemical model, boundary con¬ 
ditions of nitrate and oxygen concentrations are taken from CSIRO Atlas of Regional 
Seas (CARS 2009, http://www.cmar.csiro.au/cars) and chlorophyll a concentration from 
SeaWiES (http://oceancolor.gsfc.nasa.gov/). The simulations were performed for a 22- 
year period. The first 13 years were run with the physics only and the following 9-years 
were run with the physical/biological coupling. The coupled model reached a statistical 
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equilibrium after 4 years and model outputs were then stored every 3 days (averaged). 


Finite-size Lyapunov exponent 


The Finite-size Lyapunov exponent (FSLE), A, is a measure of the rate of divergence in 
the positions of particle pairs while separating from an initial distance (io up to a final 
distance 5f. It was developed to study non-asymptotic dispersion processed and to 
quantify dispersive behavior of particles, especially in those cases where length scales 
are easier to identify than temporal ones. It is given by the following expression: 


\ ^ 1 
A = -log^, 

r do 


where r is the time needed for the initial separation to increase from do to d/. The 
FSLE is a function of the initial and final separations, and also of the initial location 
of the particle pair xq and of the time of release Iq. Thus, the computation of A for a 
given set of initial locations and in a time interval provides an insight to the locations of 
weaker/stronger particle dispersion and its evolution with time in the domain D. In fluid 
flows, regions that exhibit substantial stretching of fluid material, hence high values of A, 
have filamental shapes and have been associated with barriers and avenues to transporP^ 
that strongly constrain the mixing of fluid with different properties, the so-called LCS 
Trajectory integration can be done from the present to the future, forward in time, 
or towards the past, backwards in time. The locations with high values of the backwards 
Lyapunov field are the structures better delimiting the distribution of transported sub- 
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stances and providing barriers to transporfi^El. Thus, the FSLE backwards field is the 
one used in this paper. 

To compute the three-dimensional FSLE field we extended a previous two-dimensional 
methocP^ to include the third dimension, by computing the time r it takes for parti¬ 
cles initially separated by 5o = (Axq + Ayl -\- to reach a final distance of 

= {Ax^i -I- Ay^i Az'jy^'^. However, in a similar application for the Benguela up- 
welling systerrpi it was observed that the displacement in the vertical .2 direction does 
not contribute significantly to the calculation of 6f and so we define a quasi-3d compu¬ 
tation of ESLE: we use the full three dimensional velocity field for particle advection 
but particles are initialized in horizontal ocean layers and the contribution Azf is not 
considered when computing 6f. 
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Figure 1: OMZ core, finite-size Lyapunov exponent (FSLE) and FSLE-O 2 gradient correlations for simulation 
year 21. a) 20 fiM isosurface of mean O 2 concentration and the mean FSLE field at 410 meters depth. Note 
the stretched vertical scale. Elat top shows the Peruvian coast, b) Zonally averaged mean (z.a.m.) ESLE and O 2 
gradient, averaged between 380 m and 600 m depth, c) Pearson correlation coefficient {R) between z.a.m. ESLE 
and O 2 gradient (solid line). Fisher 95% confidence interval on R (dashed line). 
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Figure 2: Entrainment of O 2 rich waters into the OMZ due to the motion of Lagrangian 
coherent structures. Color codifies O 2 at 410 m depth and the lines are the 0.075 day~^ 
FSLE isolines, a) 16 September; b) 7 October; c) 25 October; all of s.y. 21. Note the 
oxygen-rich tongue entering the OMZ at 80—82°W. White continuous line is the 20 yM 
mean isoline at 410 m depth (corresponding to the southern OMZ boundary). 
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Figure 3: Vertical profiles of FSLE and O 2 eddy fluxes from 200 to 600m. Mean FSLE 
(blue lines) and normal O 2 eddy flux (red lines) averaged from the coast until 88'^ W 
for each depth and during simulation year 21 at (a) Northern 20 (iM mean boundary 
and (b) Southern 20 jiM mean boundary. Positive fluxes bring O 2 into the OMZ. The 
vertical line indicates the zero value for the fluxes. 
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Backward FSLE three-dimensional fields 

Here we give additional details on the calculation of the finite-size Lyapunov expo¬ 
nent (FSLE) fields and further discuss their horizontal and vertical structure. Max¬ 
imum values of FSLE computed backwards in time can be interpreted as fronts of 
advected tracers, since they delineate the lines along which the tracers are stretched 
and folded by the fluid flow. By the same reason, they provide barriers to transport 
with little or no flow across them [1, 2]. 

We computed daily three dimensional {3d) fields of backward FSLE for simu¬ 
lation year 21. This was done by setting up a regular 3d latitude-longitude-depth 
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grid of initial conditions with spacing in the two horizontal directions of Sq = 1/27^ 
(i.e. ~ 4 km) and 30 vertical layers with variable vertical spacing. The grid covered 
the eastern tropical south Pacific (ETSP) from 88^ W to 70^ W and 18^ S to 2^ 
N. Vertically, the grid extended from 10 to 1100 m depth and was clustered in the 
upper 500 m of ocean. To obtain the FSLE field, A(x, t), at location x and time t, 
one particle is released at time t from grid node at location x and the separation 
to the particles released from neighboring nodes is monitored, r is the time needed 
for the first of these separations to reach the value Sf (that in this study was set 
at 100 km). Trajectories were integrated backward in time for 6 months with a 
Runge-Kutta 4^^ order method. If at the end of this interval the separations are 
smaller than or if the particles leave the domain or hit the shore, then the FSLE 
for the release location is set to zero. Otherwise, the FSLE is computed as 

1 S f 

A(x, t) = - log y. 

T do 

A map of instantaneous backward FSLE is shown in Fig. SI for different depths. 
In the upper layer (Fig. SI a) we observe that south of 4^ S the FSLE field is 
organized in thin filamental features with high FSLE values superimposed on a low 
FSLE background. North of 4^ S, we see a larger density of small-scale features, 
that represent the change in dynamical regime [3] as we approach the equator. Time 
scales associated with the Lagrangian dynamics at this depth vary from 4 to 10 days 
(the FSLE is roughly an inverse time scale). At larger depths (Fig. SI b, c and d) 
the features of the instantaneous FSLE field are maintained although the intensity is 
reduced from its values at the surface. The area near the equator with dense small- 
scale structure is reduced at 113 m depth and qualitatively different and shifted 
south at 410 and 592 m depth (Fig. SI b, c and d, respectively). At all depths 
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Figure SI: Map of backward FSLE for the 1** of January of simulation year 21. a) 
10 m depth, b) 113 m depth. C) 410 m depth. D) 592 m depth. Note the different 
colorbar at different depths. 
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Figure S2: Vertical map of backward FSLE for the of January of simulation year 
21, at 16.45^ S. 

shown, the high FSLE lines are more common near the coast, which in this region 
is the main source of mesoscale variability due to instability of coastal currents and 
the associated upwelling regime [4, 5]. 

The vertical structure of the instantaneous ESLE field is shown in Eig. S2 for 
the same date on a vertical section at 16.45^ S. The thin filamental structures ap¬ 
pearing in the horizontal maps have a vertical extension of about 400 m from the 
surface down. The structures oriented parallel to the zonal axis reveal this curtain¬ 
like shape, as can be observed between 76^ W and 78^ W and 100 and 300 m depth. 
Such shape has already been observed (in 3d) in similar calculations for the Benguela 
upwelling zone [6] and can be understood theoretically [7]. Dynamical studies with 
particle trajectories showed that the most intense curtains can be related to La- 
grangian eddy boundaries. 
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Figure S3: Horizontal maps of backward 
year 21. a) 10 m depth, b) 113 m depth 
the different colorbar at different depths. 
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FSLE temporally averaged for simulation 
C) 410 m depth. D) 592 m depth. Note 
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Figure S4: Time mean of the horizontally averaged depth profile of the backward 
FSLE field < A > (in day~^) for year 21. 

The temporal-mean horizontal maps of the backward FSLE fields are shown 
in Eig. S3 for four depths: 10, 113 and 410 and 592 m. One can observe the 
most persistent patterns of mixing, in particular, the northern and southern strips 
delineating the boundaries of the oxygen minimum zone (OMZ). The depth profile 
of horizontally averaged mean ESLE field for simulation year 21 (Eig. S4) shows 
the decrease in stirring with depth, with a subsurface peak at ~ 30 m depth. This 
decrease of stirring with depth, as measured by the ESLE, was also found in a 
study of the Benguela upwelling region [6]. Its origin could be simply the overall 
smaller velocities found at deeper layers, and also the decrease in the nonlinearity 
of the mesoscale eddies (as indicated for example by its ratio between rotation and 
propagation velocities) that occurs with depth [5]. 
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Residence times 


Residence times (RT) are a tool to study fluid exchange between two regions. Al¬ 
though the OMZ does not constitute a geographically determined region such as a 
basin or bay, it is still possible to deflne appropriate limits to the OMZ (as we have 
done all along this paper) so that fluid exchange between the OMZ and its exte¬ 
rior may be studied with RT distributions [8, 9]. They are a Lagrangian diagnosis 
complementary to the FSLE methods. The RT of a fluid particle is defined as the 
time it spends inside a certain region before crossing a particular boundary. Here 
we compute the RT in the OMZ as the time the particle remains with a O 2 content 
below 20 jiM. The results obtained are similar if instead we compute the exit time 
from the average OMZ core region, i.e. the volume in which the temporally averaged 
O 2 concentration is smaller than 20 fiM. Trajectories can be integrated forward in 
time (and then the computed time is properly an exit time) or backward in time (so 
that one is calculating then the time the particle has been in the region before the 
present moment). 

We show in Fig. S5 the RT distribution at 410 meters depth for the of July 
of simulation year 20. Similar results are observed at other depths. To obtain this 
distribution, particles were released from the same horizontal regular grid described 
in the previous section and the trajectories were integrated backward and forward 
for 2.5 years (the need to use such large integration times forced us to use as initial 
time the simulation year 20 instead of year 21 as in other calculations, since the 
simulation dataset contains only 22 years). Backward and forward FSLE fields were 
also computed for this date and location for comparison purposes. 

The RT are larger than 2.5 years in most of the area identified as the OMZ 
core. In contrast with this rather homogeneous distribution in the centre, at the 
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Figure S5: Residence times (RT) at 410 m depth for R* July of simulation year 
20. a) Backward RT. Superimposed white lines are 0.035 day~^ backward FSLE 
isolines. b) Forward RT. Superimposed white lines are 0.038 day~^ forward FSLE 
isolines. The magenta lines bound the region inside which the concentration of O 2 
was smaller than 20yM already at the initial time, so that outside these lines the 
exit times are zero. The red color associated to a residence time of 2.5 years is in 
fact associated also to all residence times larger than this duration. 
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OMZ boundaries the RT distribution exhibits quite complex shapes, with several 
re-entrant zones of low RT and sharp transitions to low RT. Additionally, there are 
several cases of thin regions of low RT intruding in to the high RT central zone. 
These sharp boundaries between high and low RT coincide in some areas with high 
FSLE values. This coincidence is a feature observed with the RT distributions 
of passive particles in a given geographical region [9]. Since passive particles are 
Lagrangian tracers, we conclude that in the areas of matching RT changes and 
FSLE lines, the O 2 content is conserved along particle trajectories and local changes 
in O 2 occur mostly due to advection. The long residence times found are consistent 
with the low ventilation regime accompanied with biogeochemical processes which 
reduce oxygen concentrations to hypoxic levels. 


Horizontal O2 gradients 

Here we enlarge our discussion on the relationship between FSLE and O 2 horizontal 
gradients. First we observe in Fig. S6 a) how the lines of high FSLE values determine 
instantaneous (at depth 410 m) fronts for the O 2 concentrations. These are located 
in the northern and southern boundaries of the OMZ as has been discussed all along 
the text. In Fig. S6 B) it is shown how the high-FSLE lines coincide with the largest 
horizontal gradients of O 2 . 

In Fig. S7 a) and b) we show FSLE and O 2 gradient maps, respectively, averaged 
temporally (for simulation year 21) and vertically in the mid-depth range of the OMZ 
(between 350 and 600 m). We do not discuss the relationship between the high values 
of FSLE and of O 2 gradients found very close to the coast, since that region has 
very distinct dynamics and biogeochemistry, and it is influenced by the presence of 
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Figure S 6: Lines of high FSLE values (> 0.04 day~^, black) superimposed on 
instantaneous O 2 and O 2 gradient fields on E* of March of year 21 at 410 m depth, 
a) O 2 concentration, b) O 2 horizontal gradient. 


10 
























coastal currents and water-column sediments interaction. The maps clearly show the 
high gradient zones north and south of a central basin with low average O 2 gradient 
that coincides with the OMZ. The northern band is located approximately at an 
interface between an eastward mean zonal flow relatively rich in O 2 [10, 11] and a 
southern adjacent westward mean flow. This flow configuration at middepth present 
in our simulation data shows similarities with ADCP data [12] taken in February 
2009 that shows the SICC (South Intermediate Counter Current) eastward flow and 
the adjacent SEIC (South Equatorial Intermediate Current) westward flow with 
velocities about 2.5 cmjs. The mean ESLE field for the same depth range (Eig. 
S7 a) shows that the band between 0^ and 4^ S is a zone of high ESLE. On the 
southern boundary, high O 2 gradients are located close to the equatorward edge of 
the southern subtropical gyre, along a zonal band below 16^ S. In this region, we 
also find high values of the mean ESLE but these are distributed along a wider zonal 
region, poleward from 14^ S. 

The center of the OMZ forms a basin of low O 2 gradients and also low ESLE. 
The low values of O 2 are due to low or inexistent ventilation (as quantified by large 
residence times) of this region coupled to the respiration of sinking organic matter. 
The low ESLE distribution is due to the absence of significant mesoscale activity 
at this depth. Although this region is located offshore of the Peru upwelling strip, 
and receives anticyclonic eddies formed by instability of the coastal currents [13], 
mesoscale activity as measured by eddy intensity decreases as we move offshore [5], 
and it is only intense in the two strips surrounding the OMZ. 
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Figure S7: Mean fields averaged over 300 — 600 meters depth and during year 21. 
a) FSLE. b) O 2 horizontal gradient. 

Physical and biogeochemical variability of O2 

Biogeochemical fluxes undergo a sharp decrease in variability around 300 m of 
depth, as evidenced from Figure 7 of [14]. That figure clearly indicates that the 
variability of physical fluxes becomes two orders of magnitude larger than that of 
biogeochemical fluxes below ^300 m, since the ratio is ^1% and lower. Therefore, 
O 2 can be considered as a passive tracer below 300 m. The ratio of variabilities 
(quantified as root mean square values) of physical and biogeochemical fluxes at 
12^S is shown in Fig. S8. The dataset used is the same as that described in the 
Methods section of the main text, and physical and biological fluxes are identified 
from the terms of physical and biological origin in the balance equation for oxygen 
concentration. 

Again, Fig. S8 highlights that below 300 m the variability in biogeochemical 
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Figure S8: Ratio (in %) of the variability in the biogeochemical flnxes versns the 
variability of the physical flux along a vertical section at 12'’S. The color bar applies 
to the upper panel, and the colors and values of the ratio in the lower panel cor¬ 
respond to values multiplied by 10 to highlight the vertical variability. The 45 
(white line) and 20 //M (red line) O 2 isolines are also shown. 


13 











fluxes becomes two orders of magnitude smaller than the physical ones. We note 
that over the whole OMZ, variability in physical fluxes remains in general larger 
than the variability in biogeochemical fluxes, but small biological changes in O 2 can 
still diminish correlations between physical stirring and O 2 . Below 300 m the ratio 
between fluxes drops below 1% so that one can consider that below this depth the 
physical fluxes strongly dominate over biogeochemical processes. This is consistent 
with our observations, reported in the main text, of strong correlations between 
physical stirring and O 2 distributions in the range 380-600 m, whereas this correla¬ 
tion is not so strong at upper leyers. To illustrate what is happening near the OMZ 
boundaries, we calculated the proportion of explained variance by the physical fluxes 
at different depths. The proportion of total variance (Var) in the O 2 change rate 
explained by the physical flux can be written (using that the total flux is the sum 
of the physical and the biogeochemical ones) as [15]: 


R 


2 

Phys 


1 - 


Var [Biogeochemical flux] 
Var [Total flux] 


( 1 ) 


R%hys shown in Fig. S8 at two depths (112m and 465m) corresponding to 
different zones of the OMZ. At 112 m, the biogeochemical fluxes, although having 
lower amplitude than the physical ones, can still contribute to generate local O 2 gra¬ 
dients which can influence advection terms. However at 465 m, the biogeochemical 
fluxes (within the OMZ) are almost two orders of magnitude lower than the phys¬ 
ical flux (see Fig. S8 Fig. 7 of [14]). Fig. S9 shows that physical fluxes dominate 
over most of the study region, but they are more intense near the OMZ latitudinal 
boundaries. This is more remarkable at 465m than at 112m, which indicates that 
the lower correlation between the FSLE field and mean O 2 gradient at the surface 
may arise from the relative larger contribution of biogeochemical fluxes to the rate 
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Figure S9: Explained variance of the physical flux with respect to the rate of O 2 
change in % at A) 112ni and B) 465ni. The white (red) thick line indicates the 
iso-oxygene contour at 45 fiM (20 fiM) 

of O 2 change there. 


Eddy fluxes of O 2 

We computed O 2 eddy fluxes by using the Reynolds transport theorem and the 
decomposition of the velocity and O 2 concentration fields in their mean and fluc¬ 
tuating components. Other physical fluxes such as small-scale turbulent diffusion 
turn our to be smaller. The eddy flux across each of the mean 20/iM north and 
south boundaries was computed by averaging the product of the fluctuating velocity 
component normal to the boundary in the horizontal and in the vertical directions 
and the fluctuating O 2 concentration component, thereby obtaining horizontal and 
vertical eddy fluxes. The sign of the normal was chosen so that a positive hor¬ 
izontal flux indicates transport towards the interior of the OMZ. For comparison 
purposes, mean horizontal and vertical fluxes were also computed from the mean 
velocity normal to the boundaries and the mean O 2 concentration. The integration 
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of the fluxes over the whole frontier gives the flow rate of O 2 into the OMZ through 
those boundaries. 

The profiles of mean and eddy horizontal fluxes across the 20 jiM OMZ bound¬ 
aries for depths between 200 and 600 meters were computed for simulation year 
21 and shown in Fig. SIO. The fluxes are averaged also along the horizontal ex¬ 
tension of the boundary at each depth. Blue lines are for the northern boundary 
and red for the southern one. Overall, the mean fluxes are much smaller than the 
eddy ones. In the northern boundary, eddy fluxes and mean fluxes are both positive 
up to about 250 m depth highlighting the oxygen supply role the Equatorial under¬ 
current and the primary southern subsurface counter-current play for the Pern-Chile 
under-current. At deeper locations the eddy flux is again positive, and the mean 
flow slightly negative (i.e. pointing outsize the OMZ). At the southern border the 
eddy flux is positive above 380 m, becoming negligible further down. The mean 
flux is positive above about 230 m, being negative or negligible below. The total 
mean horizontal O 2 flow rates (Table SI) at the northern and southern boundaries 
are negative whereas the eddy flow rates are positive. At both OMZ boundaries, 
mesoscale turbulence tends to transport O 2 into the OMZ at a much higher rate than 
the mean flow removes O 2 from the OMZ. At the southern boundary computed flow 
rates are higher than at the northern boundary probably due to the longer horizon¬ 
tal extension of the southern boundary. The vertical eddy and mean O 2 flow rates 
(Table S2) turn out to be about one thousand times smaller than the horizontal 
ones, so that the discussion above for the horizontal fluxes applies indeed also to 
the total fluxes. Vertical fluxes show also magnitude differences between eddy and 
mean values. 

To further investigate the correlation between FSLE and O 2 eddy flux, we singled 


16 



a 


b 




Figure SIO: Vertical profiles of horizontal O 2 fluxes at the 20fiM OMZ mean bound¬ 
aries between 200 and 600 m depth. The fluxes have been averaged along the hori¬ 
zontal extension of the boundary (from coast to 88^W) for each depth, and during 
simulation year 21. a) Eddy flux profile (these are the red curves in Fig. 3 of the 
main text), b) Mean flux profile (note the much smaller values as compared with 
the eddy fluxes.) 
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Boundary 

Eddy {jimol s ^) 

Mean [jimol s ^) 

Northern 

Southern 

6.15x10® 

1.00x10^ 

-2.01x10® 

-3.56x10® 


Table SI: Horizontal eddy and mean O 2 flow rates across the 20 //M mean boundaries 
from 200 to 600 m depth and from coast to 88^W, averaged over simulation year 21. 


Boundary 

Eddy [jimol s 

Mean {nmol s 

Northern 

1197.0 

16.5 

Southern 

-4226.4 

-269.7 


Table S2: Vertical eddy and mean O 2 flow rates across the 20 /xM mean boundaries 
from 200 to 600 m depth and from coast to 88^W, averaged over simulation year 21. 

out the contributions due to intraseasonal variability (timescales between 2-100 days, 
computed from the departure from monthly mean in the considered year 21) as 
opposed to the departures from the total annual mean used before. Indeed isolating 
the sole contribution of the O 2 eddy flux due to intraseasonal variability appeared 
interesting since FSLE might not capture stirring at the long seasonal timescales. 
In doing so, one finds an improved correlation between FSLE and O 2 horizontal 
eddy flux both at the northern and southern OMZ boundaries (linear correlation 
coefficient 0.47 versus 0.27 (north) and 0.79 versus 0.77 (south)). 


Episodic ventilation events 

We characterize the eventual entrainment of water across the OMZ boundaries by 
displaying in Fig. 11 Hovmoller plots of the time evolution of the O 2 anomaly at the 
North (panel a) and South (panel b) OMZ 20^M boundaries at 410 m depth. The 
anomaly is the actual O 2 concentration minus the mean at that location (20/uM). 
Since the overall eddy fluxes point towards the OMZ interior, positive anomalies 
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indicate oxygenated water entering into the OMZ. The slopes in the Hovmoller 
plots indicate overall westward propagation. 

O 2 anomalies are more variable for the northern (standard deviation of 4.8 jiM) 
than for the southern boundary (standard deviation of 3.1 jiM). At the Northern 
limit the O 2 front is more stable meaning that there is an available pool of O 2 rich 
waters readily available to be entrained into the OMZ core by propagating mesoscale 
structures. At the southern limit the O 2 anomaly seems to be dependent on the 
actual O 2 content of travelling eddies that occasionally approach the boundary. 

At the northern border the anomalies last longer than at the southern one and 
remain more coherent during their propagation offshore. In the South, ventilation 
events are more intermittent. However, during the year under analysis, the greatest 
episode of O 2 anomaly crossing the OMZ at this depth occurred in the southern 
boundary with peak O 2 anomaly of ^ 17 jiM around day 300 and between 600 and 
800 km from the coast (this is the episode depicted in Figure 2 of the main text). 


Cross-wavelet spectra 

As an alternative methodology to quantify correlations between O 2 concentrations 
and the stirring measure provided by FSLE we conducted a cross-wavelet analysis 
between the O 2 concentration and backward FSLE times series from our simulation 
dataset at three particular locations in order to identify the dominant time scales 
of CO- variability. The locations were: BNDNORTH_OXYCLINE, located on the 
northern OMZ boundary (84'’W,3‘’S) and 75 m depth; BNDNORTH_CORE, located 
in the OMZ core (84'’W,3'’S) at 410 m depth; and BNDSOUTH.CORE, in the 
southern boundary (84'’W,16'’S) at 410 m depth. 
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Figure Sll: Hovmoller plots of O 2 concentration anomalies at OMZ boundaries 
at 410 m depth during simulation year 21. a) Northern boundary, b) Southern 
boundary. Anomaly measured relative to 20 jiM level. Distance along boundary 
increases towards offshore, coast is at the right and the anomalies propagate offshore. 
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The climatological wavelet spectra consist in a wavelet decomposition of the 
intraseasonal anomaly time series followed by a calculation of the climatology of the 
wavelet power coefhcient at each frequency. Wavelet spectra of O 2 and backward 
FSLE signals were computed following [16] at the three geographical locations. The 
time span of the signal is 4 years, from 1st July of simulation year 18 to 30 June of 
simulation year 22. Intraseasonal anomalies were estimated as the departures from 
the monthly mean. Wavelet power was normalized following Eq. 14 of [16] in order 
to compare the magnitude of the wavelet power at different frequencies. The mother 
wavelet used was the Morlet wavelet [16]. 

The climatology of the resulting cross-wavelet power is shown in Eig. S12. The 
plots indicate where the energy is dominant, i.e. where is a significant covariance 
between O 2 and ESLE is present, as a function of calendar month. Large energy is 
found in the intraseasonal frequency band with significant peak energy around 45 
days. The marked seasonality of the intraseasonal activity indicates a link with the 
seasonal modulation of the baroclinic instability of the coastal currents (i.e. the Peru 
undercurrent) during Austral summer and in the South during Austral winter[17]. 

Analyzing the dependence of the O 2 signal variability with stirring variability, 
i.e. ESLE levels, is out of the scope of this study. This could in principle be achieved 
as in [18], namely by artificially decreasing the strength of the non-linear terms in 
the momentum equation and measuring the resulting change in the variability in 
O2. 
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Figure S12: Climatological normalized cross-wavelet spectrum of O 2 concentration 
and FSLE time series at (a) BNDNORTH.CORE, (b) BNDSOUTH.CORE and 
(c) BNDNORTH_OXYCLINE. The horizontal axis indicates the calendar month. 
Units of the colorbar are 0.8 x day~^. The white thick contour in all panels 

indicates the 95% confidence level. 
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